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We derive an exact propagation scheme for nonlinear Schrodinger equations. This scheme is en- 
tirely analogous to the propagation of linear Schrodinger equations. We accomplish this by defining a 
special operator whose algebraic properties ensure the correct propagation. As applications, we pro- 
vide a simple proof of a recent conjecture regarding higher-order integrators for the Gross-Pitaevskii 
equation, extend it to multi-component equations, and to a new class of integrators. 
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I. INTRODUCTION 

Nonlinear Schrodinger equations play a special role in 
physics. At a fundamental level, it has been suggested 
that a weak violation of the linearity of quantum me- 
chanics might allow new physics to emerge. Theoreti- 
cal studies have a long and continuing history [H, 0, [1| j 
and have resulted in several precise experimental tests 
of linearity using neutrons 0, [E @| and ions 0, [!]■ Of 
course, concrete evidence for such nonlinearity has not 
been found. At a practical level, nonlinear differential 
equations have found numerous applications in nonlinear 
optics, plasma physics, molecular dynamics in biology, 
and fluid dynamics 0. 

Most spectacularly, the dynamics of the quantum mat- 
ter wave in dilute Bose-Einstein condensates (BECs) has 
been described by the time-dependent Gross-Pitaevskii 
(GP) equation [l^: 

,f,^l^ = + W + aMr, *(r, t). 

(1) 

Among the novel phenomena successfully observed in 
BECs are four- wave mixing [ll|, solitons [H, [H, [l^ . 
and vortices [3, [l^ . Each has been accurately modelled 
by the GP equation ([T]), or its multi-component general- 
izations. 

It is natural to believe that the novel dynamics is a 
direct result of some intrinsic complexity in the time- 
evolution. In this paper, we consider the general set of 
nonlinear Schrodinger equations of the form (fi, = 1): 



(2) 



where the Hamiltonian Hq is a linear, time-independent 
Hermitian operator, the nonlinear term Kii(\E') is a real 
wave-function dependent (local) potential energy opera- 
tor, such that the normalization of 5" is preserved, and we 
have employed a simplified notation that suppresses the 



space or spin components of (see also below). We find 
that the exact propagation of Q can be accomplished 
by defining a surprisingly simple nonlinear operator Vni, 
and is given by the familiar expression 



*(t) = exp -ii(i/o + V;i) *(0). 
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(3) 

Thus, while the resulting evolution of ([2]) may be com- 
plex, the underlying rules are quite simple. 

The motivation for our work is both practical and 
fundamental. At the practical level, by understanding 
the underlying structure of nonlinear evolution equations 
such as ^ we can design more sophisticated numerical 
methods. These methods can be much more efficient the 
more we know about the formal solution. We note two 
recent examples. 

For real time propagation, it was recently conjectured 
by Javanaincn and Ruostekoski ( JR) that all higher- 
order split-operator algorithms for linear Schrodinger 
equations can be directly applied to the Gross-Pitaevskii 
(with the same order of accuracy), using only a simple 
"most-recent-update" rule to evaluate the nonlinear po- 
tential. We prove this conjecture for all orders, and for 
multi-component equations in any dimension. 

For imaginary time propagation, standard higher-order 
split-operator methods can become unstable, due to the 
presence of negative time-steps. A new class of fourth- 
order, positive time-step algorithms has recently been ex- 
tended to the GP equation in a rotating harmonic trap 
by Chin and Krotscheck (CK) [l^. Using our knowledge 
of the underlying integration properties, we show how 
to implement the positive time-step algorithms (in real 
time) for all trapping potentials. Physical results using 
these methods will be presented elsewhere. 

At the fundamental level, our work addresses the 
question: from where does the complexity of nonlinear 
Schrodinger equations arise? Our presentation shows 
that while the complexity emerges with time, it is not 
from the short-time evolution. This is done by bridging 
the gap between the functional Lie theory of partial dif- 
ferential equations known in the mathematics literature 
and the operator picture that is more common in the 
physics literature. This simple picture explicitly demon- 
strates that there is no extra computational complexity 
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to the propagation of nonlinear equations like the GP 
equation. 

The structure of this paper is the following. In Sec- 
tion II we review the propagation law for linear, time- 
independent Hamiltonians. In Section III, we present our 
main result, the propagation law ([3]) for nonlinear equa- 
tions, illustrated by the Gross-Pitaevskii equation. In 
Section IV, we prove the propagation law by introduc- 
ing the Lie operator formalism, and review the Hamil- 
tonian structure of nonlinear Schrodinger equations. In 
Section V, we study the split-operator methods for real- 
time propagation, and prove the JR conjecture, and its 
extension to time-dependent and multi-component equa- 
tions. In Section VI, we extend our results to algorithms 
with positive time steps, and conclude in Section VII. 



II. LINEAR PROPAGATION 

First, we recall that the linear Schrodinger equation 



(4) 



with time-independent Hq can be exactly propagated by 
the scheme 



(5) 



To prove (O we can either take the time derivative to 
show that ([3]) is satisfied, or we may use the Taylor series 
expression 



°0 f7l 



(6) 



n=0 



By iterating ^ we find that the n-th time derivative of 
* is 

a,"* = (7) 
and thus © can be evaluated as 

which is precisely what is meant by ([5]). 

III. NONLINEAR PROPAGATION 

For the nonlinear Schrodinger equation 

idt^ ^(Ho + Vnli^))-^ (9) 

the correct propagation is, at first sight, quite complex. 
Recall that in writing ^ we use a simplified notation 
that suppresses not only the spin and spatial coordinates 
of ^P, but also the fact that the local nonlinear potential 



Vni{^) = Vni{^, **) depends on both ^ and In this 
paper we will consider only real potentials which also 
satisfy the following reality condition: 



(10) 



This is somewhat more general than the restriction that 
the potential depend only on the local absolute value 
of i.e. T4i(\l/, = For multi-component 

equations (to be considered in Section V), where has 
components 0„ and the nonlinear potential has matrix 
elements Vnijk, this condition generalizes to 



E 



dVnLjk 



dVnlJk 



0. 



(11) 



Returning to the propagation of (O, we know that a 
Taylor expansion such as ([6]) still holds: 



(12) 



but a simple form for 9" 4' such as (O does not seem 
readily available. 

If we treat V{t) = Vni{'^{t)) as a time-dependent po- 
tential, we find i9"^' by repeated application of ([5]): 

df-^ = -i(aty)* - {Ho + V)2* 
df-^ = ~i{d^V)'^ - 2{dtV){Ho + V)^ 

-{Ho + V){dtV)^ + i{Ho + V)^^ 
df-^ = -i{dfV)'^ - 3{dfV){Ho + V)-^ - 3(9* V")^* 

+3i{dtV){Ho + F)2* - {Ho + V){d^V)-i' 

+2i{Ho + V){dtV){Ho + V)'<S/ 

+i{Ho + Vf{dtV)^! + {Ho + Vf^. 

(13) 

It appears that the number of terms in 9"^ grows with 
n, by way of the time derivatives of the potential, which 
must finally be evaluated by differentiating again. Any 
hope of achieving an expression as simple as ([7]) appears 
small. 

However, due to the deep Lie structure to the non- 
linear Schrodinger equation (to be described in the next 
section), a much simpler formulation can be given. That 
is, there exists a single nonlinear operator Vni such that 



a^vf = {~ir{Ho + KO"* 



(14) 



and whose properties we now describe. 

By comparing ^ to with n = 1, it is clear that 



Vni^ = V,a{-^)^, 



(15) 



where we have "hatted" the operator on the left to indi- 
cate that as an operator it is nonlinear in 4", while the 
un-hatted operator on the right can be considered a ^- 
dependent linear operator. Note also that does not 
uniquely specify Vni- That is, we must still define the 
higher powers of Vni and its products with Ho-, each of 
which occur in the expansion of ([1] 
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For example, a naive application of (jlSp to Hq^ would 
yield 

Vr^lHo^ = Vnl{Ho^)Ho^. (16) 

As indicated by the question mark, this result must be in 
error, for purely dimensional reasons. That is, Hq^! has 
different units than Thus, any definition should distin- 
guish between 5* and operators on 5*. Furthermore, our 



Should any Kj itself be nonlinear, its action is defined 
by using (jl7p recursively. Note that we must also specify 
the set JC. For the nonlinear Schrodinger equation, we 
require that /C = {c, Ho,Vni} and their linear combina- 
tions, where c is any complex scalar. This specification 
of K. uniquely identifies Hq = Hq x Hq, which might 
otherwise be ambiguous. 

The definition p7|) for Vni is the main result of this 
paper. Vni is a manifestly nonlinear operator, due to its 
dependence on 5" and its algebraic structure. Neverthe- 
less, it has much in common with linear operators, due 
to the following remarkable properties. First, the compo- 
sition rule is dimensionally correct, and commutes with 



Conversely, by assuming these properties we could derive 
the composition rule (fTT]) . 

Finally, note that if the function Vni satisfies 
K;(e-*^^"'(*)^') = KiK*) (which is true for the Gross- 
Pitaevskii form Vni{'i') = .g|^P), then Vni satisfies the 
power identity 

1C*=[K; (*)]"*, (22) 

in which case 

g-iAT>„,^^g-»Ay„,(*)^_ (23) 
Using the rule (|17p . the exponential propagation law 



definition should be sufficiently simple in form, yet pow- 
erful enough to generate all of the correct terms needed 
to satisfy (fT4|l . 

The correct definition for Vni is achieved in the follow- 
ing way. Let Kj, j — 1, 2, • • • n, be an element of some 
set K. of operators, to be specified. The action of Vni is 
given by the composition rule: 



(18) 



scalar multiplication: 

VnlC^ = cVnl(<f)^. (19) 

Note that c here is to be treated as an operator, and thus 
we use Vniie-'^"^) = y„,(e-*^'=^', e*^^^'*) in which 
by the condition pO)) leads to p9)) . Second, it is operator 
linear: 

VnliKi + K2)^ - VnlKi'f + VnlK2^. (20) 

Third, Vni satisfies the conjugation identity 



(21) 

I 

([3]) can be verified by simply taking its time derivative: 

idt^{t) - {Ha + exp (-it{Ho + Vni)) *(0) 

= [Ha + Ki(e-**(^°+^"''*(0))]e-"(^°+^"')*(0) 
= [Ha + Vni{-^{tm{t). 

(24) 

In deriving ()24l) . we have used the properties of the ex- 
ponential function, operator linearity (j20p and the con- 
jugation identity (PTjl . 

This derivation of the exponential propagation law, 
while correct, is only an implicit demonstration that the 
rule P?|) generates the higher derivatives An ex- 

plicit demonstration, for the Gross-Pitaevskii equation, 
is reproduced in the Appendix, where we show that 

^{t) = exp (~it{Ha + Vni)) *(0) + 0{t'>). (25) 



Vni {Ki--- Kn'i) = i" ^^^ ^"^^ [Vniie-'^'''' • ■ • g-*^"^" \I/)e-*^i • • • e-'^"^"^]^^^._;,^^o , 
where we recall that 

Vni{e-'^'^' ■■■e-'^"^"'i>) ^ Vni{e-^^'^' . . . g-'^-^-^-, e^^^^i • • • e*^"^"**). 



I 

Vnie~'^'^' ■ ■ ■ e-*-^"-^"* = Vniie-'^'^' ■ ■ ■ e-*^"-^"^)e-*^i^i • • ■ e"'^-^"*. 
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In the next section we will provide an independent argu- 
ment to show that this is correct to any order in t and 
for any nonlinear Schrodinger equation. 



IV. LIE OPERATOR PROOF 

To explicitly prove (fT4|) for general n, and therefore 
the exact exponential evolution ([3]), we use the Lie oper- 
ator formalism, which easily handles nonlinear evolution 
equations, and has found great use in the mathemati- 
cal study of differential equations [l^, [23|- In fact, all 
linear and many nonlinear Schrodinger equations can be 
written as the Hamiltonian evolution of a classical field 
[U, mi m m, HI, and t^s evolution (in an infinite- 
dimensional phase space) can in turn be described us- 
ing Lie operators [2^, [20|. These formal methods will 
justify the presentation of the previous section. While 
there have been other general studies of nonlinear wave 
equations using Lie operators [13, [H, [H, [sO], the new 
approach developed here directly shows how to convert 
the Lie differential structure into a novel algebraic struc- 
ture, a nontrivial conceptual shift that can simplify many 
calculations. 

The pair of differential equations 



28)) represents 



(26) 



(note that we consider Vni{'^*) 
be succinctly written as 



V:m^Vnim can 



where the Lie operators C = Ch 
derivative operators, defined by 



(27) 



Cv are functional 



5** 



(28) 



(30) 

with a similar expression for (j29p . The form of this ex- 
pression can be understood by introducing the Poisson 
bracket [2ll 



{AB} 



5A SB 



SB SA 



(5*(r) (51'*(r) (5*(r) (5**(r) 

(31) 

where A and B are functionals of ^'(r). This bracket 
introduces a symplectic (i.e. canonical) structure on a 
phase space where ^'(r) and ^'*(r) play the role of coor- 
dinates and momenta. The dynamical evolution of any 
phase-space function(al) is generated by a Hamiltonian 
functional, which for the Gross-Pitaevskii equation is 



d^r 



**(r)iJo*(r) + ig|vI/(r)r 



(32) 



For example, the time evolution of \E'(r) is given by 

9t*(r) = {^'(r),H} = /:*(r), (33) 

where we implicitly defined the Lie operator, i.e. C — 
{.^n} = -{n, •}. Substituting H oi §^ in the bracket 
([3T|) will produce the Lie operators jCh ([28]) and Cy ((29l) 
withVni{^)=g\^\^. 

This quick review shows how nonlinear Schrodinger 
equations can be understood as the Hamiltonian dy- 
namics of a classical field. Furthermore, linear 
Schrodinger equations correspond to Hamiltonians which 
are quadratic (or bilinear) in the field variables. We now 
show that this Hamiltonian dynamics can be represented 
by an exponential propagation law. 

First, we find that 9"^ can be elegantly written in 
terms of £ = Ch + jCy- 



(29) 



The Lie operators (|28p and ([29]) arise naturally by 
treating 5* as a classical field. Note that, in coordinates. 



(34) 



This is proven by induction. That is, we write the right- 
hand-side of ([M)) as some function F„(\E', ^t*), which we 
assume is £"5". We then compute 



= -i{Ho^ + Ki(*)*)a*F„(*, **) + *(ffo** + T4i(*)**)a*.f„(*, **) (35) 

= /:f„(*,**) 

I 



Then, using (|34p in the Taylor expansion (jl2p . we find This demonstration in fact directly implies our earlier 

^{t) = (exp(t£H + t/:y)*)^^^(o) ■ (36) 
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claim dSl), but to show this we must convert the Lie op- 
erators Ch and Cy into the original operators Hq and 
Vni- We recall a key property of Lie operators: for any 
function we have [25l. |26| 



We also define a set of Lie operators corresponding to our 
original set K, of Hermitian operators, i.e. 



d 

dm* 



(38) 



Then, since 
identity 



all Lie operators act linearly, we have the 



(37) £1 • • • CnCy-^ = 



dXi ■■■Xn 



Ai=-=A„=0 

(39) 



2 with each Ki a 



To proceed, consider the case n 
linear operator. In this case we have £"5' = ( 
In addition, Cy"^ = T^;(^')vE', and using ([37]) we compute 



^ e^i'^iT4((e^^'^^*)e^^'^^* 

= e^'^'Ki(e"*-^'-^'*)e~*-^^^^* (40) 
14;(e"*^2^2e^i^i5')e~*^2^2e^i^i5' 



Thus, the Lie operators Cj act from left-to-right, while the calculation to general n we find 
the linear operators Kj act from right-to- left. Pursuing 



I 



Substituting ([1T|) in ([55)1 and comparing with pT|) , we 
see that we should make the following identification: 

£1 • • • CnCv"^ - (-i)"Kiif„ • • • Ki^. (42) 

Finally, we observe that our assumption of linear Kj was 
not essential and can be dropped, so long as we recur- 
sively use (1171) . Using this correspondence we can now 
conclude that (|14p and ((34)) are equivalent. 

To summarize, we have shown how the Lie formalism 
directly yields a simple expression for the 5"^' ([34]) and 
exponential propagation ([36)1 . We have also converted 
the formal Lie products into a simple operator expres- 
sion by (|42p . Thus, we can associate every term in the 
Lie operator expansion of (|36p with every term in the 
expansion of ([3]). Each term is found to be identical us- 
ing the rule given by p7p . Altogether, we have shown 
how the nonlinear Schrodinger equation can be exactly 
propagated, to all orders in t, by a simple exponential. 

V. SPLIT-OPERATOR ALGORITHMS 

Recently, Javanainen and Ruostekoski (JR) have stud- 
ied the properties of higher-order split-operator ap- 
proaches to the Gross-Pitaevskii equation [l3|. A split- 



operator scheme approximates an exponential operator 
^\{A+B) ^ product of the (hopefully) simpler individ- 
ual exponentials, e^^ and e^^, in the following way: 

^X(A+B) ^ -Q gAa,AgA6,B ^ 0(A"+1), (43) 

k=l 

where the coefficients and hk and the number Nn have 
been chosen to achieve the specified order of accuracy. 
Such an approximation is known as an n-th order split- 
ting scheme. These schemes have been extensively stud- 
ied in the case that A and B are linear operators, and 
many properties of this expansion have been established. 
A review of splitting methods can be found in [3l| . 

Much less familiar is the case when either A or B 
is a nonlinear operator. It was observed in [iTj that 
the standard split-operator approach works for the one- 
dimensional Gross-Pitaevskii equation with A — —iHq 
and B = —iVni{'^) = — i.gj^'p, if the most recent update 
of the wavefunciton is used to evaluate B. For example, 
the second-order splitting (due to Strang 33]) is 

gA(A+S) ^ gAA/2gASgAA/2 ^ Q^y^y (44) 

To use this to propagate the Gross-Pitaevskii equation. 
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the correct rule is 

(45) 

where $i = exp(— iri/o/2)^'(i) is the most recent wave- 



function. JR found that this prescription will also work 
for the fourth-orde r sp litting (due to Forest and Ruth 
[13] and others [3 S [sl, [s^ ) : 



J 



,\(A+B) ^ ^XsA/2^\sB ^\(l-s)A/2^X(l-2s)B ^\(l-s)A/2^XsB ^ 0{\^) 

I 



(46) 



with s — 1/(2 — 2"'^/'^), and conjectured that this was true 
for all higher-order algorithms (note that fourth-order 
schemes for nonlinear equations were previously studied 
in [38^). 

The results of the previous section show that this 
conjecture is indeed true for all split-operator schemes, 
for (nearly) all nonlinear Schrodinger equations, and for 
multi-component wavefunctions of any dimension. That 
is, we have shown that the exact time-propagation of 
the GP equation is an exact exponential. While Vni is a 
nonlinear operator, it satisfies the important property of 
operator linearity. Thus, any split-operator scheme for 
linear operators can be directly applied to the nonlinear 
Schrodinger equation. These propagation schemes are 
particularly simple should the nonlinear operator satisfy 
the power identity ((221) • 

The explicit argument is the following. We have proven 
the exponential propagation law 



(47) 



Then, using with A = —it, A = Hq and B = Vni we 
find 



^(t 



fe=i 



To evaluate the exponentials, we use the conjugation (121 
and power ((22)) identities, such that 



where 



n 

j<k 



(50) 



By iterating (|^5)) down to fc = 1, we find that the split- 
operator scheme (|48)) becomes 



+ ^) = n e-*^"'=^°e-''^'"=^"'(*'=%(0 + 0(t"+1), 
fc=i 

(51) 

where $fc is an auxiliary wavefunction given by the 
"most-recent-update" rule: 

"J-fe = n e-'''''^"°e-'^''^^"'^'^^'><b{t). (52) 



Note that no assumptions regarding the dimensional- 
ities of Hq, Vni, and ^ have been introduced. The only 
property that is somewhat special to the GP equation 
is the power identity (|22p . Any nonlinear Schrodinger 
equation satisfying this identity can be easily integrated 
using the split-operator scheme (|?T1) - ([5^ . Previous split- 
operator approaches to the nonlinear Schrodinger equa- 
tion [39, 40, 41], have typically used only low order split- 
tings such as dm . The first important studies of higher- 
order methods include Bandrauk and Shen (ssj . who 
treat the nonlinearity as a time-dependent potential, and 
McLachlan 27] who applied Lie operator methods more 
generally. 

To highlight the generality of this approach, we apply 
it to two multi-component examples. First, consider the 
coupled nonlinear Schrodinger equations studied in [38j: 



idtc, 
idt(, 



+g2l|01p +ff22|02n 02- 



(53) 



By treating 0i and (j)2 as elements of a two-component 
wavefunction this equation takes the same form as ^ 
but where Hq and Vni are now matrices: 



Ho 



Hi 
H2 



(54) 



and 



VnlW = 



+5l2|02|^ 

521 101 



522 I 



(55) 

One can easily show that the corresponding nonlinear op- 
erator Vii satisfies the power identity ([H]), and therefore 
the integration scheme presented above can be imple- 
mented directly, since each matrix is diagonal. 

A nontrivial example is the following four-component 
equation: 



«9t0i = {Hi + Ui)c^i + 250^0304 

idt<p2 = (H2 + C/2)02 + 2.90^0304 

idt^3 = {Hi + C/3)03 + 2.9010102 

i9t04 = {H4 + C/4)04 + 250^0102 



(56) 



Equations of this form arise in four- wave mixing in BECs 
Note that the potentials Uj can depend on 0, but 
take forms similar to previous example; they are known 



Kquatii 
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as the phase-modulation terms, while the mixing terms 
have been written out explicitly. A direction matrix for- 
mulation of this equation does not obviously lead to an 
operator Vni that satisfies the power identity How- 
ever, one can split the nonlinear interactions into sepa- 
rate non-commuting terms which individually satisfy the 
identity. That is, an equivalent evolution equation reads 



VI. POSITIVE TIME-STEP ALGORITHMS 



+ Vn + Vi + V2 + V3 + V4]'^ 



(57) 



where and Vo(^') are matrices similar to ([51]) and 
and 









02'/ 


>4 



















402 

































4 




^3 






























502 














(58) 



(59) 



^^3(*) = g 



^4(*) = g 











(j) 





0401 


Vo 






















Vo 


03 01 











»1'P3 





(60) 



(61) 







Each of these matrices represents a rotation between two 
components of with a rotation angle determined by 
the other two components (which remain unchanged). 
As nonlinear operators, each satisfies the power identity 
(|22p and thus their exponential has a simple implemen- 
tation. Therefore, an integration routine for this set of 
equations can be immediately implemented, using simple 
generalizations of splitting formulae such as e.g. 

\ll{t^T) = e-*^-H"o/2g-irVo/2g-4TV'i/2g-irV'2/2 
y^^-irV3/2^-iTVi^-iTVz/2^-iTV2/2 
Xg-^ryl/2g-^ryo/2g-^r//o/2^(i) + 0{t^). 

(62) 

Finally we observe that this approach treats the non- 
linear operator Vni as a time- independent operator. For 
problems which have an explicit time-dependence (but 
not through ^f), one can use standard sequencing tech- 
niques (as indicated in [sll, see also [11]) to evaluate the 
operators. However, as observed by JR, if Vni is treated 
as a time-dependent operator, such sequencing does not 
produce expressions as simple as (l5T]) - ([5^ . 



A peculiarity of the general splitting (|43)) is that, for 
methods of third-order or higher, negative time-steps 
must appear, i.e. one of each of the coefficients {oj} and 
{hj] must be negative [3,14^14^. These negative time- 
steps have two other negative effects: they introduce po- 
tential instability in imaginary time evolution, and gen- 
erate large error terms. By introducing the new operator 
[B, [B^A\] into the product can one find positive time- 
step algorithms of third and fourth order 0, , which 
typically have nicer errors than schemes like Forest-Ruth 
[431] . One especially simple product is 

^X(A+B) ^ gAB/6gAA/2gA2B/3gAA/2gAB/6 ^ o(;^5) (g3) 



with 



B = B 



—X^BAB^A]]. 
48 L ,L , JJ 



(64) 



Recently Chin and Krotscheck (CK) [18| have shown how 
to use this factorization for imaginary-time propagation 
of Bose-Einstein condensates in a rotating harmonic trap. 
In their approach, Vni is treated as a time-dependent 
potential whose action is evaluated by a self-consistent 
iterative procedure. 

Letting A = — V^, i.e. the kinetic energy operator 
(with h = 2m = 1), and B = {Vext + Vni), and using the 
composition rule (jl7p . we find that the relevant commu- 
tator for the Gross-Pitaevskii equation, when acting on 
^f, is 



[B,[B,A]] gHV\^\^) ■ (Vl^-p) + 2.g2|*|2v2|4'|2 
+ 2g\^\^V^Vext - (VT4.t) • iS7Vext). 

(65) 

Using this expression in (j63p and (j64p , we can again con- 
struct an integrator using the "most-recent-update" rule 
just as in the previous section. Furthermore, this ap- 
proach is somewhat simpler than the self-consistent iter- 
ation method used by CK and can be applied to any 
potential, though at the cost of computing additional 
derivatives. 

There is an important distinction, however, between 
this approach and the methods developed by CK. The 
approach presented here is not directly applicable to 
imaginary-time propagation. That is, simply consider- 
ing the transformation r = it in ([2]) leads to an equation 
that does not preserve the norm of the wavefunction. 
While an exact exponential propagation law such as ([3]) 
for such an equation does exist, it requires a composi- 
tion rule slightly different from (|17[) . In addition, what is 
really implemented by CK's method is an integrator for 
the following nonlinear and nonlocal equation: 



a*(r, r) 



dr 



with 



2rn 



+ Vext{r)+M{-^,r) 



|*(r. 



jr^d^r'\^{r',T)\ 



*(i", t) 
(66) 

(67) 



8 



Development of explicit integration schemes for this op- 
erator is left for the future. 



for numerical methods but also for analytical approxi- 
mations; applications to wave-mixing in BECs will be 
presented elsewhere. 



VII. CONCLUSION 

In this paper we have presented a new approach to 
the propagation of nonlinear Schrodinger equations. By 
translating the functional Lie theory of differential equa- 
tions to a nonlinear operator, we have shown how to 
encode the propagation law as an exact exponential, 
just like linear Schrodinger equations. This formalism 
was then applied to the construction of higher-order 
split-operator methods. In particular, the "most-recent- 
update" conjecture was proven to any order, and new 
methods were proposed for multi-component equations 
arising in nonlinear pulse propagation and four- wave mix- 
ing of BECs. Finally, it was shown how these meth- 
ods could simplify the implementation of a more efficient 
fourth-order algorithm for real-time propagation in arbi- 
trary potentials (albeit not for imaginary-time). As the 
propagation of nonlinear Schrodinger equations remains 
an important topic for many physical applications, we 
believe this approach can be a productive tool not only 
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APPENDIX 

In this Appendix we show how the rule (|17p correctly 
generates the higher derivatives 9" 5* via 



(A.l) 



for n — 1 to n = A. We consider the important example of 
the Gross-Pitaevskii equation, which has Vni{'^) = 
(recall that we have suppressed the position dependence) . 
The d"'i' of (fT5|) take the complex form: 



J 



(A.2) 



2ig\^\^H^^ + 2ig^*{Ho'i' f - 4i5^'(iJo*)(i?o**) + ig^^iJo** 



(A.3) 



ig' 



-f2.gi7o(|*pi?o*) + 25i/o(**(i?o*)') - 4gi/o(*(ifo*)(-H"o**)) 
+gHo{'S^H^^*) + 2g\^\^H§^ + 6g**(i?o*)(-f^o*) 
-6.g(i7o*)^(i?o**) - 6g*(i/o**)i?o* + 65^'(i7o«')i?o** 
-g^^Hi^* + g^Hlm^^) + 2g^Hom'Hom''^)) 



V^fo(*'ifo(|*|'**)) + 6g^-^*{Ho'i')Hom^^) - 6.g2*(i/o«'*)i?o( 
-6.g2*(i/o«')i?o(|*P^'*) + 6.g2|1'|2^'2^2^* _ 4:g^Ho{\-i'\^^^Ho'S*) 

''2|^'|2i72(|^|2^) _ £,2^2^2(1^12^*) ^ Sff^i/o ( | * |*i/o^') 

-2g2|^-|27jg(^2^^^*) _^ 4.g2|*|2i/o(|5-|2iyo^) + ig^\^\2^^/*{Ho■^)^ 

-.2|^|4^2^ + 4g2^f3(i^„^*)2 _ 22g2|^r|2$(Hp^)(iJg^*) 

+g^^^Ho{'if*^Ho'i') ~ 2g'^^^Hoi\-^\^Ho-^*) ~ 252$2^o(|^'pi?o**) 
Vi^od*!'^') + 253|*|2i/Q(|$|4^) + 3g^\-^\^Hom^-9) 
-llg^\^\'^^'^{Ho^*) + 6<?3|^r|2,52^^(|^|2^*) „ c,3^'2i7o(|*|4^'*) 



,Jr|2^) 



-2.g2 



(A.4) 



To compare these with (lA.ll) . we construct the following 
non-trivial products of Vni and Hq , using the composition 
rule (HZD: 



Kii/o* = 2g\-^\^Ho^ - g-^^Ho^* 



(A.5) 
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(A.6) 
(A.7) 
(A.8) 
(A.9) 



+2g|*|2ij3^' - 65(iJo*)^(i?o**) + 6.g**(iJo*)(fi'o*) 



(A.IO) 



KiiJo'Ki* = -g2^r2ij2(|^|2^,) ^ 492*(iJo1')Fo(|*|'**) 
+ 2c/2|5r|2^2^2^* _ V^f(iJo5r*)//p(|5r|2^) 

+2g2|5-|27j2(|^|2^) _ 2g2|^|2^*(7yQ^)2 
-V|5»|2^(iJo5»)(iy(,^'*) + 4g2^f*(iJo5»)iyo(|*|25') 



(A.ll) 



+ 2.g2|^r|24»2^2^* 252^(ijQ^)iy(^(|^|2^*) 

-2.g24'(Fo**)iJo(|4'|^«') + 252^f3(iJo^*)2 
+4.g2|^r|2i7„(|^|2^(j^) _ 2.g2|4'|2iyQ(^2^^^*) 

-652|^r|2v[»(i/o^r)(iJo^r*) + 252^f*(iJo5r)ijQ([5,[2^) 



(A.12) 



V;2ij2^ = 252^-3(^^^*)2 2g2|^f|2^2^2^* _ 12^2 | 1 2 4, (//g^) (i/^^* ) 

+6g2|5'|2^*(iJo5')2 + 3,g2|5»|4iy2^ 



(A.13) 



VnlHoV^l^ = -g3^2^0(|^r|4^*) + 4g3|^-|2^2^„(|^|2^*) _ 4,g3|^f|4^2(^^^*) 

+2g^m^Hom^^) 

V^lHoVnl"^ 2g^|*|2^'2^p(|^|2^*) _ 4g3|^|4^2^^^* ^ 3g3 | ^ 1 4^^ ^ | ^ 1 2^>) 

I 



(A. 14) 
(A.15) 



Using (fOjl - jXTS)) in (|XT|) . we do indeed reproduce the 
exact expressions for 9"^' found in (|A.2[) - (|A.4p . Thus, 
we have exphcitly shown that, for the Gross-Pitaevskii 
equation, the definition for Vni (|17p and Taylor expansion 
of the exponential yields 

= cxp (-it(ifo + Vni)) *(0) + O(t^). (A.16) 
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